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Abstract 

We study the time evolution and driven motion of thin liquid films lying 
on top of chemical patterns on a substrate. Lattice-Boltzmann and molecular 
dynamics methods are used for simulations of the flow of microscopic and 
nanoscopic films, respectively. Minimization of fluid surface area is used to 
examine the corresponding equilibrium free energy landscapes. The focus 
is on motion across patterns containing diverging and converging flow junc- 
tions, with an eye towards applications to lab-on-a-chip devices. Both open 
liquid- vapor systems driven by body forces and confined liquid-liquid sys- 
tems driven by boundary motion are considered. As in earlier studies of flow 
on a linear chemical channel, we observe continuous motion of a connected 
liquid film across repeated copies of the pattern, despite the appearance of 
pearling instabilities of the interface. Provided that the strength of the driv- 
ing force and the volume of liquid are not too large, the liquid is confined to 
the chemical channels and its motion can be directed by small variations in 
the geometry of the pattern. 
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1 Introduction 



A major theme in recent technological development has been the focus on 
miniaturization and integration. The prototypical example is the evolution of 
microelectronic devices for which speed-up and cost reduction complement 
consumer appeal as the main commercial driving forces for this develop- 
ment. Analogous considerations have spawned the field of microfluidics, 
i.e., the technology of assembling complex chemical processes into a sin- 
gle miniature laboratory, as a realization of the so-called lab-on-a-chip con- 
cept Currently available devices are made of micron-scale structures but 
further miniaturization down to the nano-scale is expected 

The transition from micro- to nano-scale design and fabrication in elec- 
tronics has required new theoretical ideas in the area of mesoscopic quantum 
mechanics Pursuing the analogy, nanofluidic developments require as 
a new element the focus on new aspects of fluid mechanics. On the nano- 
scale the macroscopic description of fluids in terms of classical hydrody- 
namic equations breaks down. Boundary layers dominate, thermal fluctua- 
tions and hydrodynamic slip become relevant, the range of inter-molecular 
interactions is comparable to the system size, and the finite size of the con- 
stituent particles has to be taken into account—"—. Although this presents 
a major challenge for further miniaturization there are many ways in which 
one can make use of nano-scale specific properties of fluids, e.g., for sorting 
and sieving of biomolecules. Ion channels in biological membranes and cell 
walls are another particularly intriguing and inspiring example 

Most commercially available fluidic devices consist of closed channels 
formed by a substrate with grooves covered by a plate. Closed channels of- 
fer the advantages of preventing evaporation and of allowing for pumping by 
applying a pressure difference between inlet and outlet. However, clogging 
is a serious problem and cleaning is difficult. Open microfluidic systems, for 
which in case of a binary liquid mixture (of, e.g., water and oil) droplets and 
rivulets are confined to chemical channels formed by hydrophilic surface do- 
mains on a hydrophobic or less hydrophilic substrate, have been proposed as 
a remedy—"— . Although one cannot apply a pressure difference between in- 
let and outlet, flow in chemical channels can be induced by capillary forces, 
i.e., by wicking into channels or due to wettability gradients by elec- 
trowetting^^"^"^, Marangoni forces due to optically electrically or 
externally generated temperature gradients, (centrifugal) body forces sur- 
face acoustic waves"^^, or shear flow in a covering immiscible fluid which 
would also prevent evaporation. 

A detailed understanding of equilibrium wetting phenomena on chem- 
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ically structured substrates is the basis for designing open microfluidic de- 
vices. In particular the morphology of nonvolatile fluid droplets on chem- 
ical channels has been studied in detail and morphological transitions have 
been observed. On linear chemical channels one encounters a transition 
from bulge-shaped droplets for large volumes and large equilibrium con- 
tact angles on the channel to elongated rivulets for small volumes and small 
equilibrium contact angles on the channel On ring-shaped structures 
there is also a droplet state which covers the whole ring^^ This richness 
in morphological phases is a result of the constant volume constraint for 
non-volatile fluids. Volatile fluids, which are in chemical equilibrium with 
their vapor phase, cannot form stable droplets on homogeneous substrates. 
Accordingly, most of the morphologies on structured substrates, which are 
stable for non- volatile fluids, turn out to be unstable for volatile fluids: if the 
pressure in a fluid configuration decreases with increasing volume the drop 
will either grow without limit or evaporate since the vapor acts as a reservoir 
with fixed pressure and chemical potential. 

However, wetting experiments with volatile fluids are easier to control 
than those with non-volatile fluids. While in the latter case one must de- 
posit a tiny but well defined amount of fluid on a substrate, in the first case 
one can tune the pressure in the condensed phase by changing the chem- 
ical potential of the reservoir. In addition, transport through the vapor is 
much faster than along the surface, such that equilibration is much faster 
and not hindered by pinning of three-phase contact lines. Equilibrium wet- 
ting phenomena of volatile fluids on homogeneous substrates are rather well 
understood—!^: the fluid forms a homogeneous wetting film which grows 
in thickness as the bulk liquid- vapor coexistence line is approached from the 
gas-side. The wetting behavior of chemically patterned surfaces has been 
studied in detail ^^'^^ and the occurrence of morphological wetting transi- 
tions from very thin to thicker rivulets have been predicted for linear chem- 
ical nano-channels^^'^^. However, only recently it has become possible to 
manufacture such structures with negligible height difference between the 
channel and the surrounding substrate 

With possible applications as nanofluidic devices in mind, here we dis- 
cuss the dynamics of quasi non- volatile droplets (i.e., liquid droplets with 
a very small vapor pressure) connected by a continuous thin film, both lo- 
calized on a branched chemical channel pattern. The pattern has the shape 
of a ring which is connected periodically to copies of itself. This is the 
simplest geometry which incorporates the key features of any potential de- 
vice: one inlet, one outlet, and the splitting and reconnection of distinct fluid 
paths. This pattern has the further advantage of retaining overall periodic- 
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ity so that particles need not be added or removed from the system, which 
greatly simplifies all particle-based calculations. We do not address dewet- 
ting on chemical channels, thus putting aside the large body of research on 
dewetting of thin fluid films on homogeneous and heterogeneous surfaces 
(see, e.g, Refs. I56l462l and references therein). 

The paper is organized such that in Sect.|2]we review the lattice-Boltzmann 
(LB) method with the focus on its application to simulations of a system 
with two immiscible fluids. In Sect. [3] we discuss the shear-cell configura- 
tion used in our calculations, in which initially one fluid rests on top of a 
ring-shaped pattern as a liquid ridge, while the second fluid fills the remain- 
ing part of the simulation cell. The applied LB method does not take into 
account thermal fluctuations. Hence the results reflect ordinary continuum 
fluid mechanics and most reasonably correspond to microscale rather than 
nanoscale systems. The simulations generally show that chemical pattern- 
ing does successfully direct fluid motion, despite a hydrodynamic instabilty 
leading to the appearance of mobile non-uniform bulges or pearls of the fluid 
which prefers the chemical channels. Furthermore, a control mechanism for 
influencing fluid motion at a branching out of fluid paths can be inferred 
from the simulations. 

In Sect. |4] we apply molecular dynamics (MD)^^ as the optimal method 
for describing nanoscale systems, in order to study the same flow geometry. 
These calculations resemble earlier ones involving a linear (strip) wetting 
pattern We again find that the fluid motion is successfully directed by 
the pattern. Remarkably, the results also resemble the LB calculations, sug- 
gesting that chemical patterns provide a very robust means of flow control. 
Although the multicomponent LB calculations presented here cannot ade- 
quately handle a liquid- vapor system, there is no such limitation for MD— , 
and we present several examples of pattern-controlled motion in open sys- 
tems. When body force driven rivulets on sufficiently large patterns are 
simulated, we find a limitation for flow control: if too much liquid accumu- 
lates at a region near a junction, surface tension forces may be overcome by 
inertia and the liquid may spill off the pattern. 

Since surface tension is a dominant force in these systems and because 
the fluid behavior at the junctions of the patterns is key for controlling liquid 
motion, in Sect. [6] we study the interfacial free energy landscape of a fluid 
droplet near a junction using the Surface Evolver software package Fi- 
nally, in Sect. |7] we present conclusions and an outlook on future challenges. 
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Figure 1: The structure of the D3219-lattice. Starting from the 
central site r the velocities {ci, ...,Ci8}, allow particle transport 
to the neighboring sites located in the planes (indicated by color) 
normal to the Ci , the C3, and the C5 direction (dash-dotted vectors). 
Ci9 = (0,0,0) corresponds to particles resting at r. Upon normal- 
ization by means of the lattice constant a and the time step A^, 
the vectors Ci, C3, and C5 form a local, right handed, orthonormal 
basis so that roc = r + Ca renders the lattice sites in real space. 

2 The lattice Boltzmann method 

The LB -method is an indirect solver for the continuum fluid dynamic equa- 
tions, which mimics the underlying kinetics of the fluid particles on a lattice 
in real space with a discrete set of velocity degrees of freedom. This method 
is attractive for the current application since a number of multiphase and 
multicomponent models exist which are comparably straightforward to im- 
plement—"—. In the following, the LB model used here is described briefly 
in order to introduce the parameters which are crucial for the fluid dynamic 
description. A detailed descussion of the model can be found in the refer- 
ences cited below. 

We use a three-dimensional implementation of the multicomponent Shan- 
Chen model (SCM>^»^»22-75 ^j^j^ ^ so-called D3219-lattice^^, i.e., a three- 
dimensional cubic lattice in real space with a lattice constant a and a set of 
19 velocity degrees of freedom Ca,0fG{l,...,19}, at any lattice site with a 
position vector r (see Fig. [T]). The state of the system at a certain time / is 
given by a set of dimensionless distributions /^(r,/), which upon normal- 
ization give the fraction of particles of species ^ at a site r with the velocity 
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degree of freedom Ca, i.e., = La fa (^5 ^he total number of ^-type 
particles within the system. Resting fluid particles are taken into account by 
non-zero values of flg{r^t) with C19 = (0,0,0), whereas |ccf^i9| ^ 0. The 
particle masses do not play a physical role in the following, thus we set 
nf = m \/s, i.e., all particles of all species s have a unit mass m. 

Within our implementation of the SCM, the temporal evolution of the 
distributions fai^^t) follows the scheme 



/^(r + CaA/,/ + AO-/^(r,0 



A/ 

1 r . .... 1 

(1) 



^ 7^(r,0-C(|ca-u^P) 



where r and r + Ccf^i9A/ refer to neighboring lattice sites, i.e., off-lattice 
particle positions do not occur. The left hand side of Eq. ([TJ provides the 
convection of on the lattice along the direction of per time unit A^, 
and the right hand side effectively models a local relaxation of towards 
a dimensionless local Maxwellian velocity distribution with a relaxation 
time scale X\ The Maxwellian distribution is centered in velocity space 
around the velocity (to be defined below) and its width ct = ^JkB^ jm 
is proportional to the thermal velocity. Upon associating the internal energy 
per particle | ^5 T with the kinetic energy ^ (a/A/)^ one obtains 



ksT I a 
V m V3 A/ 

which implies T oc mk^^ (a/At)'^, so that by construction the SCM is isother- 
mal— with the temperature T fixed by the unit particle mass m, the time 
step At, and the lattice constant a. 

The zeroth moment of the distributions (r , /) leads directly to the local 
number density of species s, 

p>,0 = ^ = 3l/^(r,0, (3a) 

where n^{r,t) = La /a (1*5 renders the number of particles of species s at 
the site r. The local velocity field u^(r,/) for species s is defined by the first 
velocity moment and the so-called Shan-Chen acceleration F^, 

u'ir,t)^^{^aCafa + ^F^^). (3b) 
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The last term is needed in order to compensate for discretization artefacts 
and to obtain Navier-Stokes equations on the fluid dynamic level The cor- 
responding total quantities are given by a summation over all types: p (r, t) = 
I,p^ andu(r,0=p"^l.p'u^ 

The key feature of the SCM is the acceleration— 

{s} a 

which is often interpreted as the local effective acceleration of a particle of 
type s at the lattice site r due to the presence of particles at the neighboring 
sites r + Ca^t. The coupling strength defines attraction (negative sign 
of ^^^) or repulsion (positive sign of ^^^) between two particles of type s 
and of type s. The so-called pseudo density W^\(r,t) = 1 — ^xp[— ^^(^,0] 
determined by the local partial particle number n^{r^t). The acceleration 
enters the evolution scheme given by Eq. ([T]) via the velocity u^, which 
determines the center of the Maxwellian distribution in velocity space, 

u^(r,0=u^ + —X' 

^l^si^a^^fa + (^A? + A^)F^]. (5) 

The expressions for u^, F^, and given by Eqs. ([3bl) , (|4]), and (|5]), re- 
spectively, are constructed such that in the large scale limit one obtains fluid 
dynamic equations with a non-ideal equation of state. (In the case of a one- 
component system for a negative sign of i.e., attraction among particles 
of the same species, the pseudo density i//^ as given above leads to a van- 
der-Waals loop.) One can show that this approach is equivalent to treating 
the acceleration F^ formally as an external field. 

In the simulations, we minimalisticly mimic the kinetics of a system 
consisting of an 6>-type fluid (oil) covering a w-type fluid (water) rivulet so 
that each fluid is in contact with certain parts of an r-type substrate (rock, 
or more generally, solid). The model is minimalistic in the sense that oil 
and water particles are distinguished by a repulsive interaction only, i.e., 
> with equal relaxation time scales for both oil and water, = = 
A. The rock particles are pinned to their lattice sites, i.e., flg{r^t) = const 
and /^-^i9(r,0 = for all times respectively. Together with bounce back 
boundary conditions this effectively results in a hydrophobic solid boundary 
if/[9 >Oand ^^^>0. 

In order to induce oil-water phase separation and thus to facilitate wet- 
ting, our choice for the coupling constants is 
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The phase separation of oil and water is driven by the parameter ^, and the 
wetting behavior for a given ^ depends on the local, partial numbers densi- 
ties of oil, water, and rock particles, p^, p^, and p^, respectively. With this 
choice of parameters, phase separation is almost complete and the concen- 
tration of oil in the water phase and vice versa is negligible (smaller than 
0.2% in weight). 



2.1 The fluid dynamic level 

In the strong segregation limit, within the water or the oil phase the fluid 
dynamic equations for the majority species s E {6>,w}, i.e., for the quan- 
tities p^ and defined by Eqs. ([3ab and (l3bl) . respectively, can be ob- 
tained via a systematic Chapman-Enskog analysis of Eq. ([T])— >2ii2Zi22. With 
dt = d jdt, di = (V)/ and w- = (u^)/ for / = 1,2,3, and upon neglecting 
terms of ((u^/c^)^^^), i.e., within the regime of sufficiently small Mach 
numbers (u^/c^)^, the continuity equation is 

d,p' + di{p'u^) = , (7a) 

and the balance of momentum reads 

mp' {dt + u'jdj) u\ ^ -djp'Sij + dj ri' {dju\ + diu]) . (7b) 

In the limit of nearly complete phase separation, the density of the minority 
phase is negligible and therefore it does not contribute to the momentum 
balance. These are the compressible Navier-Stokes equations, which imply 
number conservation and conservation of the total momentum as the sum of 
all momenta of all particles (Upon (u^/c^)^ ^ the incompressible regime 
is approached asymptotically because the Chapman-Enskog analysis renders 
both Vp and V • u^o be of ((uVcr)^) .) 

For a single component and single phase system the local pressure can 
be calculated using the equation of state p ^ mc^p = ksT p of an ideal 
gas. However, for a binary mixture and for the choice of parameters given 
in Eq. ^ non-ideal contributions to the pressure come into play and one 
obtains ^ 

p^+p^+^^ y ^^VV- (8) 



mcr a -^r ^ 



With the choices given by Eq. Q and a coupling parameter ^ strong enough 
to induce nearly complete oil-water phase separation, the partial pressure 
of the minority species is negligible. The partial pressure of the majority 
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species s in the homogeneous ^-phase resembles that of an ideal gas, 

p'-p'kBT, (9) 

because a vanishing pseudo density i//^ of the minority species s renders 
the non-ideal part in Eq. ([8]) negligible. Non-ideal behavior appears only 
within the interfacial region separating the two distinct homogeneous phases 
(with majority species w and o, respectively), where the Chapman-Enskog 
analysis is not applicable. This implies that interfacial tensions and contact 
angles have to be determined from simulation results. At present there are 
no explicit formulae which allow one to calculate them directly from the 
simulation parameters. 

In a nearly completely phase separated system the shear viscosity in the 
homogeneous ^-phase reads 

ri^ri'^mp'c\ (t-^) ^t (10) 

with the dimensionless relaxation time x — Xj bd. We note that, first, rf 
depends linearly on the local number density p\ which is a property of dilute 
fluids, second, that changes upon varying the timestep Ar, and third, that 
rf does not depend on the coupling parameters (because in Eq. ^ the 
self-coupling is set to zero). 

3 Lattice Boltzmann simulations 

The chemical pattern is implemented via a lateral variation of the fixed rock 
density leading to channel-like hydrophilic surface domains with a certain 
width w embedded in a hydrophobic substrate with contact angles flphn < 90° 
and 0phob ^ 140°, respectively. (The contact angle is defined as the angle of 
the oil-water interface at a homogeneous rock surface for a macroscopicly 
large water droplet.) We have investigated two types of patterns: a symmet- 
ric pattern with xz-mirror symmetry, and an asymmetric pattern with broken 
xz-mirror symmetry (see Fig. IZja)). The simulation box is designed as a 
shear cell, which means that a no-slip boundary condition is imposed at the 
substrate, both lateral boundary conditions are periodic, and a constant flow 
velocity aligned parallel to the x-direction is imposed within the top layer of 
lattice sites in the simulation box, i.e., u^'^li'^Ol top layer ~ (^shear,0,0). 

For both the symmetric and the asymmetric pattern, we have run simula- 
tions for two different hydrophilic contact angles, 0phii ^ 30° and flphn ^ 80°, 
and for two different shear velocities, Vshear/^r = 0.12 and Vshear/^r = 0.17. 
The initial fluid geometry is given by a water rivulet covering the pattern 



9 



Figure 2: (a) We have used a simulation box with dimensions 
128^2 X 64a x 32a in x-, _y-, and z-direction and with the substrate 
in the x_y-plane as the bottom wall. The chemical pattern provides 
hydrophilic areas, marked in red, embedded in a hydrophobic sub- 
strate, marked in blue. The clipping of hydrophilic and hydropho- 
bic domains is sharp, i.e., the transition between the red and the 
blue domains occurs within one lattice spacing a. In the symmet- 
ric case all the hydrophilic sections have a uniform width w = 6a, 
and the outer radius of the ring is 20a. In the asymmetric case the 
outer radius of the front branch (the branch with smaller values of 
y) is narrowed by one a. The pattern junctions are termed the up- 
stream and the downstream junction with respect to the direction 
of the shear velocity u^'^ l^^p j^y^j. = (vshear ^ 0, 0, 0) in x-direction. 
(b) Morphology of a water rivulet on a symmetric pattern at early 
stages (100 AO for Sphn ^ 80^ 0phob - 140^ and Vshear/^r ^0.12, 
displayed as the isosurface ^ of half the density in bulk water. 
For a description of the initial configuration (t = 0) see the main 
text. The onset of the pearling instability is already visible in the 
form of small bulges above the pattern junctions. 
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Figure 3: Cuts of the density profile of the configuration shown 
in Fig. Ob), along (a) and perpendicular (b) to the xz- symmetry 
plane. The color code refers to the water density in units 
of with a bulk value (i.e., in the interior of the rivulet) of 
^ O.Ma~^. The water density drops from the bulk value to 
almost zero within a range of about 4a. The same holds for the oil 
density; however, the bulk value of the oil density (i.e., far from 
the substrate) is p"" ^ 0.74a"^, 0phii ^ 80° and 0phob - 140°. 
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with a square cross section w x w and the rest of the simulation box filled 
with oil, i.e., for the box dimensions given in Fig. |2] one has Y jw^ ^32 for 
the water volume Y . Both water and oil have homogeneous initial number 
densities, p^'^ = P = Q.l a~^, and the spatial variation of the densities be- 
tween the water and the oil phases occurs within one lattice spacing a. The 
coupling strength within the acceleration (Eq. @i) is ^ = and the 

dimensionless relaxation time is T = A/Ar = 1.0, which leads to an oil- water 
surface tension 7^ 0.07m/A^ and a shear viscosity T]^'^ = mpc\X/2 for 
both oil and water (see Eqs. (|2]) and (fTOl) for T = 1). 

With T —\ the Reynolds number and the capillary number for the water 
rivulet are Re mp^wvshear/ = (w/a) (vshear/cr) and Ca T]^ Vshear/7 = 
[{cTXfp^/{2yX^)] (vshear/cr) = [m / {^1^ J X^)] (vshear/cr), i.e., due 
to T = 1.0 and w/a 6, 

Re ^20.8^^, Ca?^0.96^^. (11) 

cj Ct 

According to our choices given above, the ratio Vshear/cr is of <^(10~^), 
so that Re is of <^(10^) and Ca is of <^(10~^). This means that we have 
laminar flow with a strong influence of the surface tension on the evolution 
of the oil-water interface. The temporal evolution of the system is followed 
by monitoring the isosurface of the water density corresponding to half its 
bulk value. 



3.1 General features 

The initially square-shaped water ridges evolve very fast (within less than 
100 At) towards a circular cross-sectional shape. Oil and water are separated 
into coexisting thermodynamic bulk phases, i.e., there is an oil bulk phase 
with a very low density of water particles (smaller than 0.2% in weight) and 
vice versa. The oil- water interface has a well-defined width of about 4a (see 
Figs.[2](Z7) and [3]). The water and oil bulk densities are somewhat above the 
initial values: ^ 0.84a~^ and p^ ^ 0.74a~^ compared to the initial val- 
ues p^ = p^ = OJa~^. This is a consequence of the Laplace pressure (recall 
Eq. ([9])) in combination with the finite volumes of the emerging interfacial 
regions. 

The water ridges are unstable with respect to pearling whereupon the 
dynamics of the pearling process for a given non-zero shear velocity v^hear 
is controlled by the values of Ca and flphu. is driven by the oil- water surface 
tension and is hindered by viscosity and the attractive fluid-substrate inter- 
actions, which means that the pearls preferentially form at the pattern junc- 
tions and the pearling is slower and less pronounced for the case flphn ^ 30° 
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compared to the case 0phii ^ 80°. Due to compressibility (recall Eq. (|7])), 
the water density in the interior decreases with a decreasing mean curvature 
of the oil- water interface, i.e., < 0.8 inside big droplets enclosing 
the vast majority of water in the system and > 0.9 inside very small 
droplets. 

3.2 Ophii ^ 80^ 

Figure m shows the time evolution for Vshear/^r = 0.12 and 0phii ^ 80°. At 
early times, the water transport due to pearling dominates over the transport 
due to shear flow. This means that whole water droplets move along the 
wetting pattern in tandem with their formation by the pearling process. The 
figure appears to show isolated droplets separated by a waterfree substrate, 
but this is an artifact of the plotting procedure. In fact, a layer of water 
with about one to two a in thickness and with approximately 1% of the wa- 
ter bulk density is permanently adsorbed at the pattern. Density fluctuation 
due to numerical round-off errors influence the fluid behavior at the sym- 
metric upstream junction, i.e., the droplets randomly choose one of the two 
branches to move downstream. At an asymmetric upstream junction, due 
to surface tension forces the wider rear branch is preferred. Bigger droplets 
sample more of the linear shear flow profile and therefore are driven faster 
than smaller ones. 

Figure [5] shows the fluid behavior for the higher shear rate Vshear/^r = 
0.17. In the symmetric case, at this rate big droplets are finally driven off the 
pattern at the upstream junction, whereas the big final droplet is still guided 
by the pattern in the asymmetric case, yet with a significant spillage onto 
the hydrophobic domain (0phob ^ 140°). Due to the enhanced shear drive, 
tails grow out of the moving droplet, which eventually break up into pearls. 
(This phenomenon might be related to the formation of Landau-Levich films 
in coating problems.) 

3.3 0phii ^ 30° 

Figures Oa) and (b) show the results for the shear rate Vshear/^r = 0.12. On 
the symmetric pattern, a single droplet gets stuck permanently at the up- 
stream junction, whereas there is no permanent hang up on the asymmetric 
pattern. The enhanced water-rock interaction leads to a water coating of the 
channel sections the droplet has run over, i.e., a water coverage far beyond 
the minimal adsorption mentioned above. 
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Figure 4: Simulation snapshots for 0phii ^ 80°, 0phob ^ 140°, 
and Vshear/<^r = 0.12 on a symmetric pattern. The pictures are 
arranged in consecutive rows with the flow of time from left to 
right. The pearling at early times, the aggregation and merger of 
droplets at the junctions, as well as the random droplet behavior 
at the upstream junction (i.e., whether turning right or left) can be 
seen. The randomness is introduced by numerical round-off errors 
in the simulation code. Note the periodic boundary conditions in 
x-direction. 
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Figure 5: Simulation snapshots for 0phii ^ 80°, 0phob ^ 140°, and 
Vshear/<^r = 0.17 on a symmetric (a) and an asymmetric pattern 
(b). In the symmetric case the droplet is driven off the pattern, 
whereas in the asymmetric case it is still guided by the wider rear 
branch of the pattern. 

For the higher shear rate Vshear/^r = 0.17 Fig. [6^0) shows that no per- 
manent hang up occurs at a symmetric upstream junction. Instead, numer- 
ical round-off errors influence the fluid behavior as already described for 
Sphii ^ 80° and Vshear/^r = 0.12. However, after the first passage the xz- 
mirror symmetry of the system remains broken by an unequal water coating, 
and hence the droplet path is predetermined for any consecutive passage by 
the film left behind the droplet. In the asymmetric case, the situation is 
qualitatively the same as the one above for Vshear/^r = 0.12. 

3.4 Flow control options 

The upstream junction where the inlet channel branches out in flow direction 
tums out to be the cmcial component of the pattem, because - provided 
that the surface tension forces are strong enough to keep the droplet on the 
pattem - the droplet behavior is determined by symmetry. Droplets get stuck 
at a symmetric junction since the mirror symmetry of the pattem naturally 
translates into a mirror symmetry of the droplet morphology, and a splitting 
of the droplet into two smaller ones is prevented by surface tension. In 
contrast, surface tension forces drive the droplets onto the wider rear branch 
of an asymmetric junction. 

This might be exploited in order to direct fluid motion and control droplet 
throughput on potential microfluidic devices: A droplet typically mnning on 
the wider rear branch of the asymmetric pattem could be occasionally forced 
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Figure 6: Simulation snapshots for 0phii ^ 30°, 0phob ^ 140°, and 
Vshear/<^r = 0.12 on a symmetric (a) and an asymmetric pattern 

(b) . While the droplet gets stuck permanently at the symmetric 
upstream junction there is no hang up in the asymmetric case. 
Due to surface tension, the droplet favors the wider rear branch. 

(c) Simulation snapshots for 0phii ^ 30° and Vshear/<^r = 0.17 on 
a symmetric pattern. The droplet at the upstream junction ran- 
domly chooses one of the branches. Due to the periodic boundary 
conditions it always takes, however, the same path in any consec- 
utive passage because the mirror symmetry of the system remains 
broken by an unequal coating of the pattern by water. The droplet 
leaves behind a significant amount of water on the channel. 
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to take the front branch by means of a time-dependent, local variation of the 
wettability near the junction generated, e.g., by electric pulses or heating. 



4 Molecular dynamics simulations 



4.1 Methods 



The MD simulations involve generic viscous fluids and a crystalline solid 
substrate constructed of atoms interacting via Lennard- Jones potentials. The 
shape of the wetting pattern on the solid is similar to the pattern used above 
for the LB simulations and the computational method is identical to that 
used in two previous papers. The case of immiscible binary liquid mixtures 
in a shear cell follows Ref. iH, while related calculations involving an open 
geometry with a liquid ridge of a one-component fluids in contact with its 
vapor is based on Ref. l39|. The pair potentials are of the Lennard- Jones 
form. 



with the atomic core size a and the potential depth e as the characteristic 
scales of length and energy, respectively. The fluid atoms of all species have 
the common mass m, so that = \fmj~k gives a characteristic time scale. 
The interaction is cut off at = 2.5<7 and shifted by a linear term so that 
the force vanishes smoothly there. The coefficient Css is used to vary the 
strength of the attractive interaction between atomic species s and s. 

In the shear cell case the interaction strength coefficients have the stan- 
dard value 1.0, except that Css = if ^ and Prefer to atoms of different liquids, 
or of the inner liquid (the "water") and of the regions outside the pattern, or 
of the outer liquid (the "oil") and of the regions of the pattern. With this, 
the two liquids are de facto immiscible and the inner liquid wets the pattern 
completely (0phii = 0°), while the regions outside the pattern are completely 
nonwetting (Sphob = 180°). 

In the open case with a one-component fluid, the liquid-liquid and the 
solid-solid interaction strength coefficients are 1 .0, and the liquid-solid coef- 
ficients are 1.0 for solid atoms within the pattern (leading to a contact angle 
Qphii = 0°) and either ("nonwetting", flphob = 180°) or 0.75 ("partial wet- 
ting", 0phob ^ 90°) for the solid atoms outside the pattern. In this open case 
the liquid consists actually of linear chains of four atoms each, joined by 
nonlinear (FENE) springs, so as to reduce the vapor pressure and to provide 
a reasonably sharp interface. 
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-6 



(12) 



17 



The atoms of the solid substrate are not free, but are tethered to lattice 
sites by a stiff linear spring with the force law f = — ^(r — Fq), where k — 
100 e/o"^ and r and Fq are the position of a wall atom and of its lattice 
position, respectively. The wall atoms have a mass = 100 m so that the 
characteristic frequency ^J~k|m^ = 1/fo of the wall atom motion is less than 
that of dimer harmonic oscillations about the LJ potential minimum, which 
is approximately 1 .6/to. 

In the shear cell case, the simulation involves 120,040 "water" atoms, 
1,799,960 "oil" atoms and 154,880 substrate atoms, in a box of dimensions 
273.6(7 X 171.0(7 in the substrate plane, with a height of 54.72(7. The of 
the wetting strip is 17.1 CJ and the outer radius of the circular arc is 68.4(7. 

The number of fluid atoms in the open flow case shown in Fig. [8] is 
119,840 and the number of substrate atoms is 84480. The box dimensions 
are identical to the those in the shear flow case. With 269040 fluid and 
345600 substrate atoms and box dimensions of 373.1(7 x 559.6(7 x 53.0(7 
the simulations leading to Figs. [9] and [TO] are roughly twice as large. 

The simulations are conducted in a box with periodic boundary con- 
ditions in all three directions, and the final liquid temperature is fixed at 
r = by a Nose-Hoover thermostat. Initially, all atoms are placed on 
the sites of a fee lattice, with the liquid on top of the wetting region having 
a rectangular cross-section; the temperature is ramped up to the above final 
value, allowing the inner liquid to assume a circular cross- sectional shape 
while remaining on top of the pattern. If left alone, the liquid forms pearls at 
the nodes of the junction pattern, due to the same surface tension instability 
as observed in the LB case. However, in the simulations discussed below a 
driving force is applied before the pearls develop. In the shear-cell case, we 
translate the atoms in the top wall at constant velocity 0.3 (7 //q, while in the 
open case a body force of magnitude 0.001 x mcr/rQ is applied to each atom 
parallel to the long x-direction of the substrate. 

4.2 Shear-cell flow 

The results of a simulation in rough correspondence to the LB calculations 
are shown in Fig. |7] The figure shows snapshots of the interface between 
two immiscible fluids in a shear cell, defined here as the surface on which 
the number density of the inner fluid is 0.4 (7~^, i.e., roughly half the bulk 
density of the inner fluid at the same temperature. The dimensions of the 
wetting pattern and of the initial volume of the inner fluid relative to the vol- 
ume of the simulation box are in similar proportions as those in the symmet- 
ric LB case in Fig. [2a). Using the definitions for the Reynolds number and 
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Figure 7: Sequence showing the evolution of a uniform hquid 
ridge (covered by a second hquid) on top of a symmetrical junc- 
tion pattern in a shear cell (see Sect. I4.1l in the main text for a full 
description of the initial configuration and the simulation param- 
eters; 0phii = 0° and 0phob = 180°). The spatial orientation of the 
sytem is the same as in Fig. Eta), i.e., the chemical pattern lies 
in the x_y-plane and the sections are aligned with the x-direction. 
The surface shown represents snapshots of the interface be- 
tween the two liquids, and the motion is driven by translating the 
top of the simulation box at a constant velocity 0.3 cj/^o along the 
positive x-direction. The frames are labeled by the time in units 
of ^0- The straight lines indicate the boundaries of the simulation 
box. 
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Figure 8: Sequence showing the evolution of a uniform, one- 
component hquid ridge on top of a symmetrical wetting pattern 
in contact with its vapor. The surface represents the liquid -vapor 
interface and motion is driven by a body force of magnitude 
0.001 mcj/^Q. The interaction coefficient cnquid.wck is 1.0 inside 
the pattern and outside; 0phii = 0° and 0phob = 180°. The other 
simulation parameters are given in Sect. 14.11 
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Figure 9: Behavior of a symmetrical junction pattern with twice 
the dimensions of Fig.O The other simulation parameters are the 
same as in Fig.[8l Sphn = 0° and 0phob = 180°. Liquid accumulates 
at the upstream node until its volume is too large to be held on the 
pattern by surface tension forces and a drop detaches and flies off. 
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Figure 10: Fluid behavior on the larger junction pattern of Fig. |9] 
but with a partially wetting exterior {cuquid.wck = 0.75); 0phii = 0'' 
and 0phob = 90°. The other simulation parameters are the same as 
in Fig. [8l In this case droplets of liquid detach from the pattern 
and slide along the substrate. 
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the Capillary number in Sec. 13.11 we have Re ^ 0.21 and Ca ^ 0.56 (com- 
pared to 2.5 < Re < 3.5 and 0.12 ^ Ca < 0.16 in the LB simulations, see also 
Sec. [5] for a discussion). The corresponding values for the viscosity of the 
inner liquid and for the liquid-liquid surface tension are rj^ ^2.l2m/{ato) 
(obtained from a separate simulation of Couette flow) and 7^ 2.50 e/<7^ 
(found from a standard simulation of a slab of liquid ^^), respectively. Fur- 
thermore, we emphasize that MD incorporates thermal fluctuations while 
(our present version of) LB does not. 

Nonetheless, the behavior of the two systems is rather similar. The 
height of the inner liquid is initially uniform, and some accumulation of wa- 
ter is visible at the nodes of the pattern before the driving force is applied at 
^ = 50^0- The fluids move from left to right (i.e., in the positive x-direction) 
in response to the top boundary motion, and at / = 1500^0 one pearl of water 
forms at the upstream node while a second one is traversing the downstream 
periodic boundary. Due to the symmetry of the wetting pattern, the bulk of 
the upstream pearl is held in place as additional water arrives through the 
periodic boundary, while small amounts advance slightly along the circular 
arcs of the pattern. As in the LB simulations, further symmetric advance 
along the arcs is impossible due to the finite volume of liquid available, and 
water continues to accumulate at the node until / = 2900 /q, after which a 
thermal fluctuation puts more material on the rear arc and the pearl moves 
along it (see the snapshot at 4100^0). Note the extreme overhang of the pearl 
- only a narrow strip along the arc is actually in contact with the solid. Also, 
as in the LB figures, the apparently water-free areas are a plotting artifact 
due to those sampling bins having less than the water density 0.4 0"~^ which 
defines the interface as an isodensity surface. At later times the single pearl 
reaches and crosses the downstream node as well as the periodic boundary 
(5800^0), halts temporarily at the upstream node again, and begins to follow 
the previous path along the rear arc (7000 ^o)- In further repeated traversals 
of the pattern at later times (not shown), the pearl configuration retains its 
shape and (approximately) its volume, and always chooses the rear arc of 
the pattern, because the thicker residue of water left behind there after the 
first passage provides a less resistive path. 

4.3 Open flows 

Similar simulations have been performed on junction patterns with only one 
one-component liquid placed on top of the wetting region, leaving the upper 
parts of the simulation box empty aside from a very dilute vapor resulting 
from evaporation of the liquid. In Fig. [8] we show the motion of the inter- 
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face of a liquid-vapor system driven by a pressure gradient. Explicitly, a 
body force is applied to each fluid atom parallel to the pattern axis, as if 
the liquid was falling vertically under gravity. However, due to numerical 
limitations the magnitude of the force had to be chosen to be much larger 
than terrestrial gravity; a practical experimental procedure corresponding to 
our simulations might use centrifugal forces on a rotating disk with a ra- 
dius of 1cm at 4 X 10"^ rpm'^^ One qualitative change (note the frame at 
lOO^o) is that the interface is rougher with much more short length scale 
variation. This may be understood intuitively as the result of molecules near 
the surface having greater freedom of movement in the absence of a viscous 
liquid "cover," and more quantitatively as a consequence of a lower surface 
tension for the liquid- vapor interface in this simulation {y — 0.49 e/o"^)^^ 
in comparison to that of the liquid-liquid interface (7= 2.50 e/<7^) consid- 
ered above. A second qualitative difference is that while in the previous 
simulation there was always a single moving pearl, more diverse morpholo- 
gies of the interface are seen here. Initially, the fluid motion is similar to 
the shear-cell case: liquid moves easily through the downstream node and 
the periodic boundary while accumulating temporarily at the upstream node 
(2050 ^o), and then randomly choosing the rear arc to traverse the junction 
(4725^0). However, on the second traversal, the pearl instead chooses the 
front arc of the junction (7425^0) because the fluid film being on the front 
arc after the first traversal through the rear arc happens to be thicker than on 
the rear arc, perhaps again due to a random fluctuation. At still later times, 
multiple pearls may be present and in different states of motion. For example 
at 11550^0, there are mobile pearls in both arcs and an accumulation at the 
upstream node of the junction. Subsequently, the bulge at the node merges 
smoothly into the two pearls which continue to move, following each other 
repeatedly through the rear arc (18750 /q)- Evidently, the larger fluctuations 
present in the liquid- vapor case lead to a less regular behavior. 

The simulations shown so far have the attractive feature that the water 
remains on top of the wetting pattern even in motion. But there are limita- 
tions to this behavior. In Fig. [9] we report on simulations of a larger system in 
which the overall size of the simulation box and the "skeleton" of the wetting 
pattern (i.e., the length of the linear sections and the radius of the centerline 
of the arcs) are doubled in size, but the width of the wetting channel and the 
initial height of the liquid are the same as before. This simulation thus in- 
volves twice the distance between nodes and twice the volume of liquid. As 
shown, in this case multiple pearls form readily, which may be understood 
from a previous analysis of liquid ridges on a linear wetting stripe which 
revealed a characteristic critical wavelength for the linear instability leading 
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to pearling behavior. Here, the longer intervals of linear or curved channels 
exceed this wavelength. A second new feature is that if the liquid is held up 
at the upstream node (due to symmetry) the additional liquid now present in 
the system may build up to a pearl too large to be held in place by surface 
tension forces until a fluctuation breaks the symmetry and an arc of the pat- 
tern is selected at random. In this situation, the pearl is simply pulled off the 
pattern by the driving forces and moves downstream in the vapor. In an al- 
ternative simulation, in which the region outside the wetting pattern features 
partial wetting (0phob ^ 90°) rather than nonwetting interactions with the liq- 
uid, as shown in Fig. [TO] streams of liquid leave the pattern and move straight 
along the substrate in the direction of the force. More generally, there is a 
competition between the wetting interactions trying to hold the liquid on the 
pattern and the applied forcing which attempts to move liquid downstream. 
If the strength of the forcing or the volume of liquid which accumulates at a 
junction is too large, the driving force pushes the liquid along the substrate 
irrespective of its wetting characteristics. 



5 Physical units 

Our simulation approach aims at reproducing classical fluid dynamic behav- 
ior via (non-fluctuating) LB and at capturing nano-scale aspects via MD. 
These two methods refer to different length scale regimes. Those regimes 
are specified by the ratio of the Reynolds and the capillary number. 

Re mpY ^^^^ 
Ca (t7^)2 ' ^ ^ 

which, for given material parameters m, p, y, and T]^, implies the value 
of the channel width w in physical units. Upon taking water under at- 
mospheric pressure and at room temperature as a reference system (i.e., 
7^ 70"^kgs"^, 7]^ ^ 10"^kg(ms)"^ and mp ^ lO^kgm"^), the ratios 
Re/Ca ^ 21.7 for LB and Re/Ca ^ 0.38 for MD (shear-cell flow) imply 
w ^ 310nm for LB and w ^ 5.4nm for MD, and hence, the spatial reso- 
lution for LB (w = 6a, see Fig. Ha)) and MD (w = 17.1 <7, see Sec. l4Jl) 
correspond to 

a^51.7nm and <7^0.32nm, (14) 

respectively. Likewise, the capillary numbers Ca 0.12 for LB (low shear 
rate) and Ca ^ 0.56 for MD imply shear velocities Vshear = Cay/T]^ ^ 
8.4m/s and Vshear ^ 39.2m/s, respectively. The intrinsic time scale for our 
systems, i.e., the typical time scale for thickness undulations of wavelength 
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w on a thin film of this thickness or of a rivulet on a channel of width w 
is given by— 

= — w = Ca . (15) 

Y ^shear 

Thus one has ^ 4.4ns for LB and ^ 77 ps for MD. The time resolutions 
of the LB and MD simulations are (see Eq. ^) At = a/{V3cT) ^ 440 ps 
and to = 0.3 a/vshear ^ 2.4ps, respectively. 



6 Free energy landscapes 

The previous analyses reveal that a key feature of motion through a branched 
pattern is the behavior of droplets at a node. In order to gain further insight 
and to elucidate the dynamics associated with this feature, we analyze a 
simplified quasi-static model which focuses on the effects of surface tension 
(see Figs.flllandO. We consider a non- volatile liquid "water" droplet near 
a node in contact both with a substrate and a fluid "oil" or vapor phase and 
perform a minimization of the interfacial free energy ^ with respect to the 
shape ^ of the liquid-fluid interface configuration for steplike variations of 
the contact angle under the constraints of a fixed volume and a fixed lateral 
position of the center of mass of the droplet (see Fig. [TJland Appendix lAb. 
For simplicity we assume the liquid to be confined to the chemical channel, 
i.e., 0phob = 180°. 

This requires to minimize the functional 

F = ^[^] + A;7Cp[^] +f|| • C|| [^]. (16) 

The minimzation constraints are encoded in the scalar expression Cp and the 
vectorial expression C|| = (Cx,Cy) with associated Lagrange multipliers Ap 
and f|| = (ix^fy), respectively. 

The constrained minimization is carried out numerically using the sur- 
face evolver package^. This is an adaptive finite element algorithm, which 
evolves a certain triangulated, initial interface configuration iteratively to- 
wards a state of minimal F via a gradient projection method with Ap, f^, and 
f^; as adjustable parameters. 

Since the droplet is confined to the chemical channel on which the equi- 
librium contact angle 9 is assumed to be constant the interfacial free energy 
term in Eq. ([T6l) has the form 

^[^] = 7 (j^dA - cos 6 • Uyj^ , (17) 
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Figure 11: The chemical pattern given by a Y-shaped (upside 
down) junction of chemical channels with a uniform width w. The 
channels are defined by a wettability jump, i.e., a steplike lateral 
variation of the contact angle 6{x^y) between hydrophilic and hy- 
drophobic values 0phii and 0phob- The channel edges are given by 
the expressions X = /(}^) = {b/2) + [y/bY-y/b] tana+w/2 
and y = g{x) = -c^l + {x/c^ tana - w[l/sina - 1/(2 tana)] 
with the parameters b — 0.2 and c = 0.2 defining the curvature 
at the comers (1) and (2), respectively. The red lines indicate 
X = ±/(}^) and y — g{x), the coordinate axes are drawn in blue, 
and a is the opening angle of the junction (here a = 45°), i.e., the 
pattern has a _yz-mirror symmetry. The origin lies in the mirror- 
symmetry plane at the potition of maximal curvature of the chan- 
nel edges, i.e.,/''(}; = 0) =0. 



26 




Figure 12: Top view (a) and bottom view (b) of the interface 
configuration ^ of a non- volatile droplet in the vicinity of the 
Y-shaped junction of chemical channels shown in Fig. [TTl ob- 
tained numerically from Eq. ([T6]) . The three-phase contact line 
^ bounding ^ is drawn in black. The droplet has a volume 
y/w^ =0.75, the opening angle is a = 45°, 0phii = 60°, and 
0phob = 180°. The lateral droplet position - indicated by the black 
cross in (b) - is specified via the center of mass projection of 
onto the x_y-plane, f y = {x^f) — (— 0.3w, — 0.8w). 

with the liquid-fluid interfacial tension y so that yw^ provides the scale for 
the free energy F . The stripe width w (see Fig. [TT]) is the only lengthscale 
within this model. Young's contact angle measures the wettability of the 
substrate, and n^; is the y-component of the interface normal pointing out- 
wards (see Appendix lAl). The substrate surface is located in the xy-plane and 
its tension with the liquid and the fluid contributes via an oriented line inte- 
gral along the boundary of i.e., the three-phase contact line ^ (Fig. [12]). 
The chemical pattem is given by a steplike lateral variation of the contact 
angle (see Figs. [TTl and[T2l). 

Non- volatility is ensured upon constraining the liquid volume to a fixed 
value y. This constraint is implemented via the Lagrange multiplier so 
that 



where JA is the oriented surface element, and is the z-component of the 
outward normal of the interface. The Lagrange multiplier A;? can be inter- 
preted as the Laplace pressure of the droplet. 

The constraint C|| is used in order to fix the center of mass projection of 




(18) 
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Figure 13: The topography (a) and the corresponding contour 
plot (b) of the free energy F{x^f) as a function of the center of 
mass coordinates x and y in the x_y-plane for the system shown in 
Fig. [121 The blue cross within the contour plot indicates the origin 
of the coordinate system (see Fig. [TT]) and the black lines corre- 
spond to the positions of the channel boundaries (see Fig. [TTI) . 
There is a minimum of F at (x = 0,_y — 0.4w). 



28 



-1 -0.5 0.5 1 

X I w 



-1 -0.5_ 0.5 1 

y I w 



Figure 14: Cuts of the free energy landscapes for V /w^ = 0.75, 
a — 45°, 0phob = 180°, and various hydrophilic contact angles 
Qphii ^ [60°, 110°], (a) along x through the global minimum and 
(b) along the jf- symmetry plane at x = 0. Due to the symmetry of 
the system the global minimum always occurs at x = 0. In (a) and 
(b) the free energy is shown relative to the free energy of a droplet 
on a homogeneous hydrophobic substrate Fsub = F{x ^ j) 
and on a homogeneous straight chemical channel Fchan = ^(1-^1 < 
^ oo), respectively. Therefore the values for the depths of the 
minima in (a) and (b) differ. 
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the droplet onto the substrate plane at a certain position ry = {x^y), i.e., 




see Appendix El The vectorial Lagrange multiplier = (fx,fy) can be inter- 
preted as a spatially constant, lateral force field, which suitably balances the 
effect of surface tension forces. In other words, the configurational space 
{^} for the minimization of the free energy functional ^[^], which has 
been restricted to a subspace of interface shapes with a certain volume V 
by the contraint Cp given in Eq. (fTSl) . is restricted further to the subspace of 
shapes with a certain center of mass projection f||. Naturally, this method 
only works for droplets with a finite support, i.e., for Qphii > 0°. 

With this, the free energy F of a droplet on the pattern can be obtained as 
a function of x and 3; with V and the wettability contrast 0phob — ^phii as input 
parameters (see Fig.[T3]for the topography and a corresponding contour plot 
of F{x^y)). Naturally, the chemical channels translate into energetic valleys 
with a global minimum at (x = 0,3; ^ — 0.4w) near the intersection point of 
the channel midlines at 3; = ^ (l — (sin45°)~^) ^ — 0.2w. Upon decreasing 
the wettability contrast the topography of the free energy landscape becomes 
less structured: the depth of the global minimum relative to the free energy 
of a droplet far from the chemical channels as well as relative to the free 
energy of a droplet on a straight channel decreases (see Figs. [Mta) and (b), 
respectively). However, the position of the global minimum does not shift 
(see Fig. [14] (b)). The y-position of the minimum certainly depends on the 
opening angle a. For a = 60°, i.e., for a threefold rotational symmetry of 
the junction, it is right at the intersection point of the channel midlines. 

The occurrence of the global free energy minimum at the pattern junc- 
tion implies a certain qualitative behavior of the droplet: An unconstrained 
droplet - given by f y =0 within Eq. ([T6l) - would be sucked in and trapped at 
the junction. Equivalently, a trapped droplet has to overcome a free energy 
barrier in order to escape from the junction. In this resect it is worthwhile 
to compare this with the LB-results in Figs. |4l|6l which show preferential 
droplet formation on the junction during pearling, as well as the hang-up 
and merger of driven droplets there. In particular Fig. 0a) shows a perma- 
nent hang-up for 0phii ^ 30°. 

The contour plot in Fig. \T3\h) shows, that the consequences of the free 
energy minimum, i.e., leading to a possible droplet hang-up, can be avoided 
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by guiding the droplet via an additional external, lateral force in such a way 
that its center of mass projection f y moves along an open contour line. 

7 Summary, conclusion, and outlook 

We have analyzed theoretically the forced motion of liquid rivulets through 
branched chemical patterns on a substrate (Figs. |2] and |3}. The purpose 
of this study is to demonstrate the feasibility of this set-up for controlled 
fluid transport, which might serve as a component of future "lab-on-a-chip" 
devices. Our calculations employ a lattice Boltzmann method (Fig. [T]) for 
the isothermal description of microscale systems and molecular dynamics 
for their nanoscale counterparts. 

The microscale simulations (see Figs.|4]-[6]) show that within reasonable 
limitations for the flow rate and the liquid volume, the wetting pattern is ca- 
pable to direct liquid motion along desired paths. Familiar surface tension 
instabilities lead to pearl formation at the nodes of the pattern (see in partic- 
ular Fig. |4]), but as in simpler flow configurations— the pearls are mobile 
and move along the pattern, and therefore need not to obstruct the motion 
of the liquid. Although the nanoscale systems are particularly sensitive to 
thermal fluctuations, the behavior of the liquid in this domain is nonetheless 
fairly similar (see Figs.rzllTOl). 

Concerning the possibility to control these systems the challenge is to 
direct the liquid preferentially along one of multiple paths leaving a node. 
Our microscale simulations indicate that small variations in the width of 
the wetting channels can suffice for this purpose (see Figs. |5jb) and Ob)). 
The liquid prefers to move along a wider channel where viscous resistance 
is lower. Other mechanisms for directing flow, such as imposed electric 
fields acting on an ionic liquid, bursts of air directed towards one side of a 
node, smooth wettability gradients on the surface, etc., might work as well 
or better. We leave these possibilities for future studies. 

For a qualitative understanding of the behavior of driven pearls at the 
nodes we have employed a quasi-static model based on the action of surface 
tensions. This yields the free energy landscape with a state of minimal free 
energy for an isolated droplet sitting at a node (see Figs. [lTI-[T4l). 

Likewise, an additional consideration for future work is the behavior of 
suspensions or complex liquids, e.g., lyotropic or thermotropic liquid crys- 
tals, flowing along chemical patterns. 
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A Free energy 



Using the local version of Young's law, 7cos0(x,3;) = Yrf{x^y) — Yrwi^^y) 
with the indices r, w, and / for "rock", "water", and "fluid" (i.e., "oil" or 
vapor), respectively, within the capillary model for a chemically heteroge- 
neous substrate domain with a laterally varying macroscopic contact angle 
9{x^y) the free energy ^ can be expressed as 



— = / dA — [ dA-h^cosOix^y) 

7 7^ Jy 



(20) 



with the liquid- fluid interfacial tension 7 as a scaling paramter. Here we as- 
sume that the substrate surface and thus also the "water"-substrate interface 
y lie in the xy-plane and that y is oriented such that its normal vector is 
= (0,0,1). 

The integral over y can be converted into an integral over the three- 
phase contact line if one finds a vector potential for the vector field cos 9 (x, y) . 
A simple (and for the present geometry convenient) choice of the corre- 
sponding gauge is adopted by 

n, cose(x,y) = V x [hyC(x,y)] , (21) 

with n^y = (0, 1 , 0) and C(x, 3;) = /^^ dx' cos (x' , 3;) . (Note that V x [n^ C(x, 3;)] 
is independent of xq.) Applying Stokes' theorem yields 



^ I dA - i d\-nyC{x,y), (22) 

where the three-phase contact line ^ is the oriented boundary of y as well 
as of (see Fig.[T2l). 

If a chemical pattern with a steplike variation of the local equilibrium 
contact angle 0(x,y) confines the drop to a chemically homogeneous part of 
the substrate Eq. (l22l) reduces to Eq. ([TTI) . This holds even if a portion of jSf 
coincides with an interval of the boundary of the chemical channel. 



B Center of mass constraint 

In order to constrain the lateral position of the center of mass Jy dV r of 
a droplet of volume V and of homoegeneous density to (x^y) one adds the 
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expression 




with the Lagrange multiplier f|| = (fxify) to the free energy functional. The 
volume integrals in Eq. (l23l) can be converted into surface integrals using 
GauB' theorem. To this end one expresses each of the scalar integrands as 
the divergence of a suitable vector field. For droplets residing on a planar 
substrate which lies in the xy-plane a convenient choice is 

x^V-C^x^n,) and 3; = V-(i/n,). (24) 

Since the area element dS = dSUz of y is orthogonal to ftj^ and the con- 
tribution of the surface integral over the liquid-substrate interface y to the 
surface integral over the total surface of V vanishes and the term • C|| in 
Eq. ([T6l) reduces to the expression given in Eq. ([T9l) . 
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